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Abstract 

The  goal  of  the  project  was  to  develop  a  new  computational  code  to  be  used  at  the  Arnold  Engineering 
Development  Center,  (AEDC),  Tullahoma,  Tennessee,  for  analysis  of  plasmas  characterized  by  high  ionization 
fractions  and  extremely  low  number  densities.  Such  conditions  arise  in  two  particular  test  facilities  operated  at 
AEDC:  (1)  the  Decade  Radiation  Test  Facility  (DRTF),  used  to  test  the  effects  of  strong  radiation  sources  on 
materials  and  spacecraft  components;  and  (2)  the  12V  Space  Chamber,  used  to  test  electric  propulsion  (EP) 
thrusters  used  to  control  spacecraft.  In  order  to  have  accurate  predication  of  the  quasi-neutral,  kinetic  flows  in 
both  facilities,  the  research  was  focused  on  developing  a  general,  unsteady,  electro-static,  nonequilibrium  gas 
and  plasma  simulation  code  that  combined  the  capabilities  of  the  direct  simulation  Monte  Carlo  (DSMC) 
method  and  the  Particle-In-Cell  (PIC)  method.  Efforts  were  also  concentrated  on  applying  the  code  to  study  the 
flows  in  the  Decade  and  12V  facilities. 

Introduction 

The  objective  of  this  project  was  to  develop  a  new  computational  code  to  analyze  the  flows  of  plasmas 
characterized  by  high  ionization  fractions  and  extremely  low  number  densities,  for  application  to  the  Nuclear 
Weapons  Effects  (NWE)  Decade  facility  and  to  the  12V  Electric  Propulsion  (EP)  plumes  facility  located  at  the 
Arnold  Engineering  Development  Center  (AEDC).  Physically  accurate  numerical  simulations  are  expected  to 
play  a  significant  role  in  optimizing  the  performance  of  the  facilities  as  well  as  in  playing  an  integral  role  in  the 
testing  capabilities  provided  by  AEDC. 

Typical  flow  conditions  in  Decade  and  12  V  give  a  Debye  length  and  a  mean  free  path  that  indicate  that  the 
flow  is  in  the  quasi-neutral,  kinetic  flow  regime.  Traditional  continuum-based  computational  codes  are 
unsuitable  for  very  low  density,  or  rarefied,  gas  and  plasmas.  Particle  methods,  such  as  the  direct  simulation 
Monte  Carlo  method  (DSMC)  [1],  are  becoming  mature  for  modeling  rarefied  flows  of  neutral  gases,  but  do  not 
account  for  charged  particle  interactions.  Particle-in-Cell  (PIC)  methods  [2]  account  for  field  interaction  effects 
of  charged  particles,  but  do  not  fully  resolve  particle  energy  and  momentum  transfer  effects.  Codes  combining 
the  capabilities  of  DSMC  and  PIC  must  be  utilized  to  account  for  all  the  afore-mentioned  effects  in  NWE  and 
EP  testing.  This  research  effort  will  build  upon  existing  DSMC  and  PIC  techniques  to  develop  a  combined 
PIC/DSMC  code  directly  applicable  to  AEDC's  test  facility  needs:  a  general,  unsteady,  3D,  electro-static 
nonequilibrium  gas  and  plasma  simulation  code  that  is  readily  usable  by  engineers  at  AEDC  for  analysis  of  the 
Decade  and  12V  facilities. 

Accomplishments 

Starting  with  an  existing,  general,  3D  DSMC  code  called  MONACO,  we  added  many  new  capabilities  required 
to  accurately  model  the  AEDC  facility  flows,  performed  solid  code  validations,  and  studied  typical  flows 
experienced  in  the  Decade  and  12V  facilities. 

Code  Development 

We  have  added  many  features  to  MONACO  that  expands  the  capability  of  the  code  and  improves  the  interface 
to  users.  Motivated  by  the  Decade  unsteady  flows,  we  added  the  following  features: 

•  The  ability  to  routinely  run  unsteady  simulations. 

•  The  ability  to  accept  multiple  flow  boundary  conditions. 

•  The  ability  to  set  multiple  initial  flow  conditions. 

•  The  subcell  technique  for  collision  pair  selection  [  1  ]. 
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•  The  near-partner  technique  for  collision  pair  selection  [3]. 

•  The  ability  to  perform  dynamic  domain  decomposition  for  the  unsteady  flows  executed  on  parallel 
computers  by  integrating  Metis  package  [4], 

•  The  relaxation  technique  to  estimate  macroscopic  flow  properties  [5], 

•  Feedback  to  users  on  simulation  parameters. 

•  A  utility  to  generate  particle  weight  of  cells. 

Many  of  these  features  improve  the  generality  or  efficiency  of  the  code.  For  instance,  the  subcell  selection 
technique  allows  the  cell  size  of  the  grid  to  be  larger  than  the  mean  free  path  of  molecules,  thus  avoids  a  new 
simulation  using  a  finer  grid.  The  parallel  domain  decomposition  technique  enables  the  workloads  to  be  equally 
distributed  on  all  processors,  which  can  improve  the  parallel  efficiency  by  several  times.  The  code  has  been 
widely  validated  and  the  results  agree  very  well  with  published  data.  The  validation  ranges  from  one¬ 
dimensional  to  three-dimensional  flows,  and  includes  subsonic  and  supersonic  flows.  In  addition,  work  this 
year  focused  on  the  development  of  a  new  condensation  model  that  is  consistent  with  the  DSMC  technique  [6]. 
This  modeling  was  found  to  be  required  due  to  the  very  low  temperatures  experienced  in  DECADE,  and  is 
described  in  detail  below. 

The  12V  flows  involve  plasma  jet  expansions  from  an  electric  propulsion  device  into  the  vacuum  chamber. 
Thus  far,  for  these  flows,  we  have  added  the  following  capabilities  to  MONACO: 

•  The  ability  to  simulate  electro-static  plasma  using  the  Particle-In-Cell  (PIC)  method.  The  ions  are  treated 
as  particles  and  the  electrons  as  a  fluid. 

•  The  general  unstructured  grid  implementation  for  the  PIC  method. 

•  The  important  collision  mechanism  of  charge  exchange. 

•  Particle  weighting  by  species. 

•  New  diagnostics  (e.g.,  ion  energy  distribution). 

More  technical  details  are  provided  in  the  later  portions  of  this  report. 

Transitions 

Updated  version  of  MONACO  software  delivered  to  Ken  Tatum,  Aerospace  Testing  Alliance,  Arnold  AFB. 
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TECHNICAL  DETAILS 


1.  Decade  Flow  Investigations 

In  recent  years  there  has  been  much  interest  in  the  behavior  of  clusters.  The  formation  and  growth  of  clusters 
through  gas  nucleation  and  condensation  is  believed  essential  to  many  phenomena.  For  instance,  vapor-liquid 
nucleation  processes  play  an  important  role  in  the  formation  of  atmospheric  aerosols.1  The  formation  of  soot 
begins  with  formation  of  clusters.2  Gas  condensation  can  also  occur  within  thruster  plumes,  and  clusters  may 
cause  contamination  of  the  rocket  surfaces.3  The  formation  of  clusters  will  affect  the  supersonic  flow  in 
expansion  nozzles  due  to  the  heat  released  during  the  phase  change  process.4  Cluster  formation  and  growth  is 
also  involved  in  many  material  processes,  including  chemical  vapor  deposition,5  dry  etching,6  and  nanoscale 
material  manufacturing.7  Understanding  the  mechanism  of  the  gas  nucleation  (cluster  formation)  and 
condensation  (cluster  growth)  will  therefore  help  improve  fabricating  techniques  and  the  quality  of  the  resultant 
products. 

Studies  on  gas  nucleation  and  condensation  are  generally  based  on  two  approaches:  continuum  approach  and 
kinetic  approach.  The  continuum  approach  is  usually  called  classical  nucleation  theory8,9  that  is  derived  from 
thermodynamics.  It  is  assumed  that  homogeneous  nucleation  starts  when  the  free  energy  loss  from  the  transition 
of  gas  molecules  into  the  liquid  phase  can  compensate  for  the  energy  increase  resulting  from  the  surface  tension 
of  a  cluster.  Different  expressions  have  been  derived  for  the  nucleation  rate,  and  some  of  them  can  differ  by  a 
factor  of  1017  for  the  magnitude  of  the  nucleation  rate.9  Furthermore,  the  surface  tension  of  clusters  is  still 
unclear  yet.  If  the  value  of  the  bulk  material  is  used  for  the  surface  tension  of  clusters,  an  advanced  expression 
(Lothe  and  Pound10)  for  the  nucleation  rate  predicts  a  much  worse  nucleation  rate  than  a  simple  expression 
(Becker  and  Doring8)  as  compared  with  measurement  data.  Therefore  there  are  still  many  uncertainties  in  the 
classical  nucleation  theory.  On  the  other  hand,  many  efforts  have  been  made  recently  to  develop  kinetic 
approaches.  One  is  called  Smoluchowski’s  approach  where  nucleation  is  viewed  as  a  process  of  chemical 
aggregation."  The  key  issue  is  to  determine  the  rate  constants  that  are  still  unresolved  yet.  A  fundamental 
approach  for  studying  nucleation  and  condensation  is  the  molecular  dynamics  method  that  simulates  molecular 
movement  and  interactions  directly.12  This  approach  is  supposed  to  expose  the  detailed  nucleation  and 
condensation  processes,  and  has  been  applied  to  study  gas  nucleation  and  condensation.8,13'16  However,  studies 
using  molecular  dynamics  are  limited  to  systems  with  few  clusters  or/and  clusters  having  small  size  because  of 
the  large  numerical  cost  for  simulations.  It  is  impractical  to  simulate  meso-scale  systems  such  as  supersonic 
flow  in  an  expansion  nozzle  based  on  the  current  computational  capability. 

A  more  efficient  kinetic  approach,  the  direct  simulation  Monte  Carlo  (DSMC)  method,17  has  been 
introduced  into  the  field  of  gas  nucleation  and  condensation.  The  DSMC  method  is  less  numerically  expensive 
than  molecular  dynamics  but  still  kinetically  accurate,  and  gas  nucleation  and  condensation  models  can  be 
easily  implemented  in  the  DSMC  method.  Therefore,  the  DSMC  method  is  becoming  a  popular  approach  for 
model  testing  and  flow  investigation.11,18'22  For  instance,  Hettema  and  McFeaters"  used  the  DSMC  method  to 
implement  the  Smoluchowski  approach;  Zhong  et  al.21  employed  the  classical  nucleation  theory  in  the  DSMC 
method  studying  supersonic  plumes.  Most  of  these  applications,  however,  have  neglected  the  cluster  effects  on 
gas  flows,  which  then  affects  the  cluster  modeling.  There  are  two  significant  influences  of  clustering  on  gas 
flows:  the  removal  of  a  portion  of  the  vapor  phase  and  the  “heating”  of  the  remainder  to  absorb  the  energy 
extracted  from  the  condensed  phase.  In  the  usual  case,  the  ratio  hfg/CpT  {hfg  being  the  enthalpy  change  during 

condensation,  c  being  the  specific  heat,  and  T  is  the  local  temperature)  is  greater  than  unity  and  this  means 

that  “heating”  has  a  larger  effect  on  the  stream  properties  than  the  vapor  removal.4  These  effects  were  partially 
demonstrated  in  the  plume  studies  by  Perrell  et  al.  although  a  continuum-based  approach  was  used.23  We  also 
showed  the  clustering  effects  on  gas  flow  in  our  previous  paper  when  studying  gas  expansion  in  a  supersonic 
nozzle  by  including  dimers  in  the  simulation.22  However,  there  is  no  systematic  nucleation  and  condensation 
model  developed  for  the  DSMC  method  yet. 

The  main  goal  of  this  study  is  to  develop  a  general  nucleation  and  condensation  model  to  be  implemented  in 
the  DSMC  method.  The  model  itself  is  not  able  to  predict  the  value  of  parameters  for  simulations,  but  is 
supposed  to  utilize  new  results  from  molecular  dynamics  simulations.  The  rest  of  the  paper  is  organized  as 
follows.  Brief  description  of  homogeneous  nucleation  and  condensation  theory  is  given  in  section  II.  Then 
microscopic  modeling  that  is  suitable  for  the  DMSC  method  is  discussed  in  section  III.  Several  parametric 
studies  are  given  in  section  IV.  Finally,  the  paper  ends  with  some  concluding  remarks. 


Homogeneous  Nucleation  and  Condensation  Theory 

A  gas  is  regarded  as  undersaturated,  saturated,  or  supersaturated  if  the  pressure  of  the  gas  ( P  )  is  less  than,  equal 
to,  or  larger  than  the  equilibrium  saturation  value  (Ps)  corresponding  to  the  local  temperature.  The 

undersaturated  and  saturated  states  are  thermodynamically  stable,  whereas  the  supersaturated  state  is  unstable. 
The  supersaturated  gas  will  nucleate  and  condense  in  the  presence  of  impurities  or  walls.  Otherwise,  the  gas 
sustains  a  very  high  degree  of  saturation  (S  =  P/PS)  until  homogeneous  nucleation  occurs. 


There  are  generally  two  different  views  on  the  homogeneous  nucleation  processes.  The  classical  nucleation 
theory  (CNT)  indicates  that  nucleation  begins  after  the  “energy  barrier”  is  overcome;  whereas  kinetic 
approaches,  such  as  Smoluchowski’s  approach,  assume  that  nucleation  starts  from  formation  of  dimers.  In  the 
classical  nucleation  theory,  it  is  assumed  that  clusters  grow  or  shrink  via  the  attachment  or  loss  of  a  single 
molecule.  Making  this  approximation  leads  to  a  set  of  coupled  rate  equations  for  the  number  density  of  clusters 
of  different  size.  Further  with  the  equilibrium  state  assumption,  the  nucleation  rate  ( J)  can  be  expressed  as 
follows: 


j_cP<$*r>)  (  A g'\ 

JmPt  P(  kT  J 

where  r.  is  the  critical  cluster  radius,  and  AG*  is  the  free  energy  of  formation  of  a  critical-sized  cluster.  The 
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other  variables  are:  Boltzmann  constant  k ,  gas  number  density  n ,  molecular  mass  m  ,  and  constant  c .  The 
critical  size  is  defined  such  that  the  free  energy  of  formation  reaches  the  maximum  value  at  this  size,  and 
clusters  larger  than  this  size  are  stable  whereas  smaller  clusters  are  thermodynamically  unstable.  The  CNT  relies 
on  a  macroscopic  approximation  for  evaluation  of  the  free  energy  of  clusters.  This  of  course  makes  no  sense  for 
small  clusters  of  a  few  molecules.  There  are  also  arguments  on  the  number  of  degrees  of  freedom  to  be  included 
in  the  evaluation  of  free  energy.9  In  general,  the  nucleation  rate  can  be  expressed  using  a  nucleation  constant 
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where  o  is  the  surface  tension  of  clusters,  and  P/ ,  pv  are  the  mass  densities  of  liquid  and  vapor,  respectively. 
The  critical  cluster  radius  is  given  as 
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Determining  the  value  of  a  is  another  challenge.  There  is  no  genera!  agreement  on  how  the  surface  tension 
relates  to  the  radius  of  curvature.4  The  surface  tension  is  also  temperature-dependent.  Although  there  are  several 
uncertainties  in  the  CNT,  measurement  data  can  actually  follow  the  CNT  prediction  (expression  (2))  for  most 
cases.8  However,  when  the  vapor  temperature  is  low,  the  critical  size  predicted  by  CNT  may  not  be  reasonable. 
For  instance,  for  condensation  of  water  vapor  in  air  at  low  temperatures  T  <  210K  ,  the  size  of  the  critical 
cluster  determined  by  CNT  is  less  than  the  water  molecule  itself.24  Some  even  suggested  not  to  use  the  CNT  for 
modeling  the  onset  of  clustering  in  supersonic  beams.25 

Kinetic  approaches,  however,  do  not  use  the  concept  of  “equilibrium  state”  and  there  is  no  critical  cluster 
size.  Clusters  are  formed  when  two  molecules  approach  each  other  and  form  a  bound  state.  It  is  generally 
assumed  that  the  formation  of  a  bound  state  is  feasible  only  when  a  third  particle  interacts  with  one  of  the  two 
particles  and  carries  away  the  kinetic  energy.  Molecular  dynamics  simulation  by  Zeifman  et  al.  confirmed  that 
condensation  starts  from  dimer  formation  in  triple  collisions  of  monomers.26  The  nucleation  rate,  however,  has 
not  been  generally  determined.  Because  clusters  can  dissociate,  the  cluster  distribution  can  diverge  as  the 
volume  containing  the  clusters  increases.  Treatment  of  the  volume  dependence  of  the  free  energy  has  been  a 
source  of  confusion  and  controversy  in  the  development  of  a  general  molecular  theory  of  nucleation.  Other 
kinetic  approaches,  however,  cannot  predict  the  nucleation  rate.  For  instance,  the  Smoluchowski  approach  treats 
cluster  formation  as  a  chemical  reaction,  but  the  value  of  rate  constants  is  unavailable  from  the  approach  itself. 
Another  difference  is  that  a  third  particle  is  not  included  for  the  cluster  formation  reaction  in  the  Smoluchowski 
approach.  Overall,  there  is  no  general  expression  for  the  nucleation  rate  in  kinetic  approaches. 

There  are  fewer  discrepancies  about  the  growth  and  depletion  of  clusters.  Both  continuum  and  kinetic 
approaches  assume  that  cluster  condensation  and  evaporation  proceeds  via  one  monomer  at  a  time.  Other 


(3) 


processes  such  as  re-arrangement  and  fragmentation  processes  are  sometimes  included.  When  a  monomer  has  a 
collision  with  a  cluster,  a  sticking  coefficient  ( q  )  is  used  to  indicate  the  possibility  for  a  condensation  event. 

Then  the  condensation  mass  flux  can  be  expressed  as  qsp/j2jtRT-  A  collision  between  a  monomer  and  a  cluster 

can  also  activate  the  cluster  to  evaporate.  The  evaporation  mass  flux  is  derived  such  that  the  evaporation  mass 
flux  is  equal  to  the  condensation  mass  flux  in  equilibrium,  and  is  expressed  as  ^  Pdf^2jiRTd  .  where  Td  is  the 

cluster  temperature  and  p  is  the  hypothetical  ambient  pressure  which  would  be  necessary  to  keep  the  cluster  in 

equilibrium. 


Microscopic  Modeling 

Gas  nucleation  and  condensation  are  essentially  kinetic  processes.  Detailed  description  of  these  processes, 
however,  is  numerically  too  expensive;  and  it  is  unnecessary  for  many  practical  applications.  A  statistical 
description  is  then  more  useful  for  numerical  simulations.  In  this  section,  we  will  discuss  some  properties  of 
clusters,  and  develop  microscopic  models  that  can  be  easily  implemented  in  the  direct  simulation  Monte  Carlo 
method  for  gas  nucleation  and  condensation. 

Some  Properties  of  Clusters 

Theoretically,  clusters  bridge  the  gap  between  microscopic  gas  and  macroscopic  materials.  The  behavior  of 
small  clusters  (dimer,  trimer,  etc.)  is  governed  by  atomic  and  molecular  mechanics  whereas  the  behavior  of 
large-sized  clusters  is  governed  by  the  macroscopic  properties  of  the  material.  Therefore,  as  clusters  form  and 
grow,  their  behavior  develops  from  the  molecular  to  the  bulk  behavior. 

Clusters  are  treated  as  particles  so  that  microscopic  behavior  of  clusters  can  be  well  represented.  For 
simplicity,  only  a  monatomic  gas  is  discussed  in  this  study.  If  we  denote  a  j-mer  as  a  cluster  having  j  atoms, 
then  the  mass  of  a  j-mer  is  jm  .  Its  translational  energy  (  Ete  )  is  jmVj  /2 ,  where  Vj  is  the  translational  speed  of 
the  j-mer.  The  average  of  the  translational  energy  is  then  3kTj/2  from  gas  kinetic  theory,  where  Tj  is  the 
temperature  of  the  j-mer.  The  internal  energy  (E,e)  of  a  j-mer  is  assumed  to  have  an  average  of  3{j -\)kTi  [2  so 
that  the  total  of  the  translational  energy  and  internal  energy  of  a  j-mer  is  equal  to  the  corresponding  total  of  j 
atoms.  Unlike  gas  atoms,  clusters  have  potential  energy  ( Erz ),  which  can  be  expressed  as  -  jhjf,  {fj ,  j  J.  Here, 
hfgMj’j)  *s  t'le  latent  heat  of  evaporation  of  the  cluster.  Clearly,  h Ifj ,  j  J  should  be  small  for  small  clusters 

because  there  are  relatively  few  bonds  for  each  atom;  and  it  should  approach  the  latent  heat  of  evaporation  of 
liquid  (hj-g(Tj))  when  the  cluster  size  increases.  The  latent  heat  of  evaporation  is  temperature-dependent,  and 

its  relation  can  be  derived  using  the  first  law  of  thermodynamics: 

Ete  +  E,c  +En+Q  =  Etlo  +  £(£0  +  Em  (4) 

|*r + 1(/  -  i>r  -  jhk  (T)+jCpk(T„  -  r).  hr,  +|(/-i>r0  -jhj t,)  (5) 

^)-[c„-fjk(r0-Thhfi(T0)  (6) 

where  Q  is  the  heat  required  to  heat  the  cluster  from  T  to  T0  .  The  dependence  of  the  latent  heat  on  cluster  size 
is  still  unclear.  If  we  refer  to  the  often  used  binding  energy  expression  Eb(j)=  aj  +  asj^} (flv  and  as  are 

constants  corresponding  to  the  volume  and  surface  terms),19  the  latent  heat  of  evaporation  of  clusters  is  assumed 
as  follows: 

En(T,j)-  -jh/s(T$-CJ-'/>  +CJ-')  (7) 

where  the  last  term  is  added  to  increase  options  for  parametric  studies.  If  clusters  are  assumed  as  spheres,  the 
radius  is  therefore  based  on  the  volume.  In  fact,  the  radius  of  clusters  is  better  represented  by 

rj-Ap  +  B  (8) 

as  predicted  by  molecular  dynamics  simulation,27  where  A  and  B  are  constants.  Molecular  dynamics 
simulations  also  show  that  the  cluster  radius  depends  on  the  relative  collision  velocity  as  in  the  variable  hard 
sphere  molecular  model:17 


(9) 


where  v  is  the  temperature  index  of  the  viscosity,  which  also  depends  on  the  cluster  size.27 


Nucleation 


As  discussed  in  Section  II,  there  is  no  general  kinetic  expression  for  gas  nucleation.  The  expression  from  the 
classical  nucleation  theory  is  then  used  to  derive  a  microscopic  expression  for  the  nucleation  rate.  A 
multiplicative  factor  may  be  used  in  the  expression,  and  parametric  studies  may  be  performed  to  study  effects  of 
this  factor.  Here,  we  assume  that  the  formation  of  dimers  is  the  onset  of  nucleation  as  in  general  kinetic 
approaches.  For  simplicity,  the  third  particle  is  also  omitted  and  the  metastable  collision  complex  is  assumed 
stabilized  upon  the  next  collision.  This  assumption  is  acceptable  since  we  do  not  really  know  what  the  value  of 
the  nucleation  rate  is.  The  nucleation  rate  is  then  converted  to  a  nucleation  probability  of  binary  collisions. 
Namely,  a  nucleation  event  occurs  according  to  the  probability  of  nucleation,  which  is  derived  as  follows: 


C  JSma/ji  (  4nr?o  \ 


where  aT  is  the  collision  cross  section. 

After  a  nucleation  event  occurs,  the  status  of  the  formed  dimer  can  be  evaluated  as  follows: 


mV,  +  mV,.  =2  mV2 

(ID 

lmy’  +imK|?  -X-2mVt  +  |  kT2  +  Ef[(t„2) 

(12) 

Vi-\(v,+Vy)>hT2-^V,'-En(T2a) 

(13) 

Condensation  and  Evaporation 


When  a  monomer  collides  with  a  cluster,  there  are  three  possible  collision  types:  condensation,  evaporation, 
and  regular  reflection.  We  combine  the  condensation  and  evaporation  probabilities  into  one  probability,  and  use 
the  sign  of  the  combined  probability  to  indicate  a  possible  nucleation  or  condensation  event.  The  combined 
condensation-  evaporation  probability  is  derived  as: 


Pc 


Here  a  positive  value  of  pc  indicates  condensation,  and  a  negative  value  refers  to  evaporation.  If  neither 


(14) 


condensation  nor  evaporation  occurs,  the  monomer-cluster  collision  is  a  scattering  collision  where  full  thermal 
accommodation  is  typically  assumed.  The  status  of  particles  after  collisions  are  derived  as  follows: 

Condensation 
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Evaporation 
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Note  that  the  post-collision  status  is  determined  only  for  a  condensation  event.  The  status  for  an  evaporation 
or  reflection  event  is  not  determined.  A  general  procedure  for  this  undetermined  status  is  to  use  the  Borgnakke- 
Larsen  model,17  which  is  explained  in  the  next  sub-section. 


Larsen-Borgnakke  Model 

The  Larsen-Borgnakke  model  is  widely  used  in  the  DSMC  method  to  redistribute  total  collision  energy 
between  translational  and  internal  energy  modes  by  sampling  energy  from  two  equilibrium  distributions.  Using 
the  equal-  partition  principle,  the  distribution  of  energy  Ea  from  a  total  of  Ea  +  Eb  can  be  expressed  as  follows: 


/(*)=/! 


E„  +  Eh 


rfe.+g»)r. 


(23) 


where  ^  and  are  the  numbers  of  degrees  of  freedom  for  Ea  and  Eb ,  respectively.  The  instant  value  of 
energy  Ea  is  sampled  using  the  acceptance-rejection  method.17  Namely,  a  value  of  Ea  is  chosen  randomly 
between  zero  and  Ea  +  Eb ,  then  this  value  is  used  to  calculate  its  probability  ratio  using  Eq.  (24)  and  is 
compared  with  a  random 
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fraction  R j  that  is  generated  from  a  uniform  distribution  between  0  and  1 .  This  value  of  Ea  is  accepted  if  the 


probability  ratio  is  greater  than  R j ,  but  a  new  value  is  chosen  and  the  process  is  repeated  if  the  ratio  is  less  than 


R  j .  However,  when  ^  and  are  not  of  the  same  order  of  magnitude,  many  attempts  have  to  be  made  to  find  a 
value  of  Ea.  A  remedy  is  to  avoid  finding  x  itself  by  converting  x  to  a  new  variable  y  whose  value  is  close 
to  0.5  by  using  a  transformation  y  -  x^a  .  The  value  of  a  can  be  chosen  as  .  Then  the  distribution  for 

y  is 


and  the  probability  ratio  is 
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To  compare  the  numerical  performance  of  the  standard  and  modified  Larsen-Borgnakke  models,  the  average 
number  of  attempts  required  during  the  acceptance-rejection  procedure  is  plotted  in  Fig.  1  (a)  where  the 
percentage  of  one 
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Figure  I.  Numerical  performance  of  the  Larsen-Borgnakke  model,  (a)  comparison  of  number  of 

attempts  during  the  acceptance-rejection  procedure,  (b)  instant  distribution  of  the  energy  percentage  using  the 
standard  Larsen-Borgnakke  method. 

energy  in  the  total  is  Ea/(Ea  +  Eb).  Clearly,  the  modified  Larsen-Borgnakke  method  improves  greatly  the 
numerical  efficiency  when  |  is  much  smaller  than  which  is  usually  the  case  for  clusters  having  many 

atoms.  Numerical  tests  also  show  that  the  standard  Larsen-Borgnakke  method  can  accept  unusual  values  with  a 
relatively  large  probability.  For  example,  Fig  1(b)  shows  the  instant  energy  percentage  (Ea/(Ea  +Eb)) 
distribution  whose  theoretical  average  value  is  2e"4.  We  find  that  there  are  4  values  larger  than  0.1  among 
10,000  values,  which  is  statistically  too  much.  These  undesired  values,  however,  are  avoided  when  the  modified 
Larsen-Borgnakke  method  is  employed  (the  plot  is  otherwise  similar  to  Fig.  1(b),  and  is  therefore  not  shown). 


Parametric  Studies 

The  developed  kinetic  model  is  implemented  in  the  DSMC  research  code  “MONACO”.28  The  clusters  are 
modeled  as  particles  in  a  similar  way  to  monomers.  Clusters  have  translational  velocity;  and  their  internal  and 
potential  energy  is  represented  by  the  cluster  temperature.  The  sub-relaxation  technique29  is  employed  to 
evaluate  the  macroscopic  properties  (gas  temperature,  pressure,  et  al.)  that  are  explicitly  involved  in  the 
developed  microscopic  model.  We  use  argon  gas  to  illustrate  the  nucleation  and  condensation  processes. 

As  discussed  in  Section  III,  there  are  many  uncertainties  about  the  nucleation  and  condensation  processes. 

For  simplicity,  the  clusters  are  modeled  as  hard  spheres,  which  means  that  the  cluster  radius  is  calculated  as 
y'^r) .  The  surface  tension  of  clusters  is  approximated  as  0.0344(1  -  T/Tc)  N/m  following  Hale,30  where  the 
critical  temperature  (Tc)  is  about  150.85K.  The  saturated  vapor  pressure  Ps  is  used  as  Pd  for  the  evaporation 
calculation.  The  saturated  vapor  pressure  is  written  using  a  logarithmic-exponential  curve-fit  expression  as 
follows: 

\nPa  =-396.465  +  225.279  In T-  41. 7722  In2  T  +  2.631  A  In3  T  (27) 

when  the  temperature  is  in  the  range  between  20K  and  1 50.85K. 

Two  numerical  examples  are  employed  to  study  effects  of  several  parameters  on  gas  nucleation  and 
condensation.  One  is  a  box  flow.  The  argon  gas  represented  by  90,000  particles  is  simulated  in  a  square  box 
where  the  box  surface  is  assumed  specular.  The  initial  number  density  and  temperature  are  lxl026/m3  and  55K, 
respectively.  Other  parameters  are  set  as  following  unless  specified  otherwise:  c  =  1 0* >  C0  =  2.0,  and 

C  „  1.175.  The  second  example  involves  supersonic  flow  in  an  expansion  nozzle. 


Nucleation  Rate 


There  is  not  much  information  about  the  nucleation  rate  at  low  temperature.  It  is  therefore  very  important  to 
study  effects  of  this  rate  on  flow  properties.  Figure  2  shows  some  plots  using  three  different  values  for  the  rate 
constant.  In  this  study,  evaporation  of  a  dimer  to  two  monomers  is  disabled  because  fluctuation  of  flow 
properties  cannot  nucleate  a  dimer  after  a  dimer  is  accidentally  evaporated.  These  plots  show  that  the  number  of 
clusters  increases  linearly  at  the  early  time  and  then  quickly  reaches  a  constant  as  the  degree  of  supersaturation 
decreases  because  the  monomers  are  converted  into  clusters  and  gas  temperature  is  increased.  The  average 
cluster  size  is  almost  independent  of  the  nucleation  rate  at  very  early  time,  but  then  increases  in  a  behavior 
dominated  by  the  condensation  process  (smaller  nucleation  rate  corresponds  to  larger  clusters).  Both  gas  and 
cluster  temperatures  increase  as  more  and  more  gas  is  converted  into  clusters,  and  the  gas  temperature  lags 
behind  the  cluster  temperature.  It  takes  more  than  a  tenth  of  a  microsecond  to  reach  the  steady  state.  In  general, 
a  larger  nucleation  rate  means  that  there  will  be  more  clusters  but  having  a  smaller  average  size,  and  the  gas  will 
be  heated  quicker  but  there  is  less  effect  on  the  final  flow  temperature  which  is  basically  determined  by  the 
vapor  pressure. 


Sticking  Coefficient 

The  sticking  coefficient  is  usually  pre-assumed  for  condensation  simulations.  Recent  studies  of  Zhong  et 
al.27  showed  that  the  sticking  coefficient  depended  on  the  cluster  size.  For  the  time  being,  we  only  consider  a 
fixed  sticking  coefficient  for  all  clusters.  Figure  3  shows  comparisons  of  several  flow  properties  using  three 
different  values  for  the  sticking  coefficient.  Clearly,  the  sticking  coefficient  has  no  effect  on  the  number  of 
clusters  at  early  time,  but  a  larger  sticking  coefficient  means  a  large  cluster  size  and  later  limits  the  number  of 
clusters  at  steady  state.  Because  the  sticking  coefficient  determines  the  condensation  rate,  it  also  affects  the 
temperature  of  both  gas  and  clusters.  A  larger  sticking  coefficient  will  heat  the  gas  more  quickly.  If  the  sticking 
coefficient  is  very  small,  then  gas-cluster  collisions  will  have  a  larger  possibility  for  a  reflection  collision;  so  the 
gas  temperature  can  closely  follow  the  cluster  temperature. 
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Figure  2.  Effects  of  nucleation  rate  on  flow  properties,  (a)  number  of  clusters ,  (b)  average  cluster 

size,  (c)  gas  temperature,  (d)  cluster  temperature. 
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Effects  of  sticking  coefficient  on  flow  properties,  (a)  number  of  clusters,  (b)  average  cluster 


size,  (c)  gas  temperature,  (d)  cluster  temperature. 
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Figure  4.  Effects  of  potential  energy  on  flow  properties,  (a)  number  of  clusters,  (b)  average  cluster 

size,  (c)  gas  temperature,  (d)  cluster  temperature. 


Potential  Energy 

The  value  of  potential  energy  is  also  unclear  for  small  clusters.  We  therefore  test  effects  of  potential  energy 
on  flow  properties.  Figure  4  shows  some  results  using  several  different  values  for  the  potential  energy.  Maybe 
because  of  the  small  sticking  coefficient  (0.1),  there  are  no  large  effects  of  the  potential  energy  tested  on  the 
flow  properties.  Further  studies  are  required  to  draw  solid  conclusions  for  effects  of  the  potential  energy. 

Supersonic  Flow  in  Expansion  Nozzle 

It  is  not  unusual  that  gas  nucleation  and  condensation  occurs  in  supersonic  flow  in  expansion  nozzles.4,21'22 
We  use  a  simple  example  to  illustrate  the  gas  nucleation  and  condensation  phenomena  in  supersonic  expansion 
flows.  The  expansion  nozzle  is  two  dimensional,  and  the  surface  is  assumed  specular.  The  inflow  argon  gas  has 
a  uniform  velocity  of  400m/s,  a  temperature  of  60K,  and  a  number  density  of  1024/m3.  Values  of  several 


parameters  used  in  the  developed  kinetic  model  are  set  as  follows:  c„  -to28,  q,  -1.0,  Ca  =2.0,  and  c  =1.175. 
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Figure  5.  Supersonic  flow  in  an  expansion  nozzle,  (a)  number  of  clusters,  (b)  average  cluster  size,  (c) 

number  density  of  clusters  (1/m3),  (d)  mass  fraction  of  clusters,  (e)  gas  velocity  (m/s),  (f)  cluster  velocity  (m/s), 
(g)  gas  temperature  (K),  (h)  cluster  temperature  (K). 

The  flow  patterns  of  this  problem  are  illustrated  in  Fig.  5.  Figure  5(a)  shows  the  number  of  clusters  formed 
in  the  simulated  domain  (the  average  particle  number  in  each  cell  is  about  500).  Clearly,  the  clusters  are  formed 
in  a  relatively  small  region  (Fig.  5(c)),  and  the  cluster  size  increases  as  the  clusters  move  downstream  in  the 
flow  (Fig.  5(b)).  So  the  mass  fraction  of  clusters  also  increases  downstream  the  flow  (Fig.  5(d)).  It  is  found  that 
the  velocity  of  clusters  (Fig.  5(f))  is  very  close  to  the  gas  velocity  (Fig.  5(e)),  whereas  the  cluster  temperature 
(Fig.  5(h))  is  larger  than  the  gas  temperature  (Fig.  5(g))  because  there  are  few  reflection  collisions  between  gas 
molecules  and  clusters.  If  we  decrease  the  sticking  coefficient,  then  few  gas  molecules  will  condense  onto  each 
cluster  (Fig.  6(b)).  However,  there  are  more  clusters  formed  in  the  flow  (Fig.  6(a)).  The  temperature  of  both  gas 
and  cluster  decreases,  and  the  difference  between  them  also  decreases.  If  we  decrease  the  nucleation  rate 
(c„=  1022  ),  then  clusters  are  formed  further  downstream  of  the  nozzle  (Fig.  7(a))  where  the  degree  of 


supersaturation  is  larger.  The  average  size  of  dusters  increases,  and  the  temperature  of  both  gas  and  clusters 
decreases  as  expected. 
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Figure  6.  Supersonic  flow  in  an  expansion  nozzle  with  a  smaller  sticking  coefficient  (0.1).  (a) 

number  of  clusters,  (b)  average  cluster  size,  (c)  gas  temperature  (K),  (d)  cluster  temperature(K). 
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Figure  7.  Supersonic  flow  in  an  expansion  nozzle  with  a  smaller  nucleation  rate,  (a)  number  of 

clusters,  (b)  average  cluster  size,  (c)  gas  temperature  (K),  (d)  cluster  temperature  (K). 


Discussion  and  Conclusion 

Gas  nucleation  and  condensation  is  very  common  in  both  nature  and  industry.  There  have  been  many  studies 
on  this  topic.  The  classical  nucleation  theory  is  able  to  explain  some  physics  of  the  nucleation  processes,  and 


molecular  dynamics  simulation  is  supposed  to  expose  the  detailed  mechanism  of  gas  nucleation  and 
condensation.  However,  due  to  the  physical  inaccuracy  of  the  classical  nucleation  theory  for  small  clusters  and 
the  extremely  expensive  numerical  cost  for  molecular  dynamics  simulations,  there  is  no  general  tool  for 
studying  gas  flows  having  clusters. 

The  direct  simulation  Monte  Carlo  method  is  a  less  numerically  expensive  but  still  kinetically  accurate 
method.  We  develop  a  general  nucleation  and  condensation  model  and  implemented  it  in  the  DSMC  method, 
and  thus  detailed  flow  simulation  becomes  possible  for  flows  involving  clusters.  In  this  model,  clusters  are 
modeled  as  particles,  but  both  internal  energy  and  potential  energy  are  considered.  Gas  nucleation  and 
condensation  are  modeled  via  particle  collisions.  The  probability  for  possible  nucleation  or  condensation, 
however,  cannot  be  determined  from  the  model  itself.  Since  there  are  many  uncertainties  in  the  physics  of 
nucleation  and  condensation,  the  DSMC  method  is  very  useful  to  investigate  effects  of  parameters  on  general 
flow  properties. 

Our  parametric  studies  showed  that  the  proposed  model  was  able  to  test  effects  of  different  parameters.  The 
nucleation  rate  affected  not  only  the  number  of  clusters  and  cluster  size  but  also  the  properties  of  the  remaining 
gas.  Investigation  also  showed  that  the  sticking  coefficient  was  very  important  for  both  cluster  and  gas 
properties.  The  application  of  the  proposed  model  to  a  model  example,  supersonic  flow  in  an  expansion  nozzle, 
demonstrated  the  capability  of  the  proposed  model:  it  is  able  to  simulate  complicated  flows  involving  clusters 
and  can  predict  detailed  flow  properties  of  both  gas  and  clusters. 
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2. 12V  Flow  Investigations 

HALL  thrusters  represent  an  efficient  form  of  plasma  electric  propulsion  for  spacecraft.  They  offer  a  high 
specific  impulse  that  is  well  suited  for  satellite  station-keeping,  repositioning,  and  orbit  transfer.  There  are 
concerns,  however,  about  the  plumes.  For  instance,  the  plumes  may  contaminate  spacecraft  surfaces  and 
interfere  with  satellite  communications.  Such  effects  need  to  be  understood  during  the  development  of  thrusters 
and  their  integration  onto  spacecraft.  Successful  integration  of  plasma  thrusters  onto  spacecraft  involves  a 
mixture  of  analysis  and  ground-based  experiments  conducted  in  vacuum  chambers.  A  key  aspect  of  the  vacuum 
chamber  experiments  is  the  desire  to  maintain  as  low  a  back  pressure  as  possible.  Elevated  back  pressures  can 
lead  to  augmentation  of  thrust  and  distortion  of  the  plasma  plume,  thereby  complicating  the  process  of 
integration  assessment.  The  12V  vacuum  chamber  located  at  the  Arnold  Engineering  Development  Center 
(AEDC)  is  a  large  facility  with  a  total  pumping  rate  of  about  3-5x1 06  litres/sec  on  xenon.  For  thruster  mass  flow 
rates  of  10-20  mg/sec,  it  is  able  to  maintain  a  back  pressure  on  the  order  of  10'6  torr.  Further  details  of  the 
facility  and  its  operation  can  be  found  in  Ref.  1 .  In  addition  to  providing  a  physical  test  capability,  AEDC  has  a 
desire  to  provide  computational  analysis  support  for  customers  interested  in  using  their  facilities.  The  focus  of 
the  present  work  is  on  the  development  of  an  analysis  too!  that  can  be  applied  at  AEDC  to  model  the  operation 
of  different  plasma  thrusters  in  the  12V  vacuum  chamber.  The  present  study  is  limited  to  investigation  of  Hall 
thruster  plumes. 

In  general,  the  near  plume  of  a  Hall  thruster  consists  of  neutrals,  highly  energetic  ions,  and  electrons.  The 
particles  have  collisions  due  to  their  thermal  velocities,  and  some  collisions  between  neutrals  and  ions  lead  to 
charge-exchange  interactions,  which  produce  slow  ions  and  highly  energetic  neutrals.  Furthermore,  the  ions  are 
also  affected  by  the  self-consistent  electric  fields.  In  addition,  a  background  gas  is  always  involved  in  the 
ground  tests  of  thrusters.  Therefore,  the  behavior  of  thruster  plumes  is  very  complicated.  An  efficient  approach 
for  understanding  these  processes  is  to  use  computer  modeling2  because  the  physics  at  different  levels  can  be 
included  for  the  plume.  For  instance,  the  direct  simulation  Monte  Carlo  (DSMC)  method3  can  be  used  to 
capture  the  collision  dynamics,  and  the  particle-in-cell  (PIC)  method4  can  be  applied  to  include  the  electric  field 
effects.  In  addition,  computer  simulations  can  identify  the  relative  importance  of  the  physics  involved  in  the 
plume. 

In  this  paper,  particle  simulation  of  a  Hall  thruster  plume  in  12V  is  investigated  using  a  hybrid  DSMC/PIC 
code.  The  rest  of  the  paper  is  organized  as  follows.  Section  II  introduces  the  numerical  method.  Section  III 
describes  the  plume  simulation.  Section  IV  investigates  effects  of  different  physics  and  compares  simulation 
results  and  measurement  data.  Finally,  some  concluding  remarks  are  given  in  Section  V. 

Numerical  Method 

Hall  thrusters  primarily  use  xenon  as  propellant.  The  plume  is  comprised  of  beam  ions  with  velocities  on  the 
order  of  20  km/s,  low  energy  charge-exchange  ions,  neutral  atoms  from  the  thruster,  electrons,  and  the 
background  gas  of  the  experimental  facility.  The  interactions  of  these  species  as  well  as  the  influence  of  the 


electric  fields  are  the  important  modeling  issues.  Computational  analysis  of  Hall  thruster  plumes  is  regularly 
performed  using  a  hybrid  particle-fluid  formulation.  The  direct  simulation  Monte  Carlo  (DSMC)  method 
models  the  collisions  of  the  heavy  particles  (ions  and  neutrals).  The  particle-in-cell  (PIC)  method  models  the 
transport  of  ions  in  electric  fields. 

Collision  Dynamics 

The  DSMC  method  uses  particles  to  simulate  collision  effects  in  rarefied  gas  flows.  The  particles  represent 
real  ions  and  neutrals,  and  are  grouped  in  cells  whose  sizes  are  less  than  a  mean  free  path.  Pairs  of  these 
particles  are  selected  at  random  and  a  collision  probability  is  evaluated  that  is  proportional  to  the  product  of  the 
relative  velocity  and  collision  cross  section  for  each  pair.  The  probability  is  then  compared  with  a  random 
number  between  zero  and  one  to  determine  if  that  collision  occurs.  If  so,  some  form  of  collision  dynamics  is 
performed  to  alter  the  properties  of  the  colliding  particles. 

There  are  two  types  of  collisions  that  are  important  in  the  Hall  thruster  plume:  elastic  (momentum  exchange) 
and  charge  exchange.  Elastic  collisions  involve  only  exchange  of  momentum  between  the  participating 
particles.  For  the  systems  of  interest  here,  this  may  involve  neutral-neutral  or  neutral-ion  collisions.  For  neutral- 
neutral  collisions,  the  variable  hard  sphere  collision  model3  is  employed.  The  collision  cross  section  for  xenon 
is: 


oEL(Xe,Xe)= 


2.12x10  2 

- ; - m 

g  “ 


where  g  is  the  relative  velocity  and  a>  =  0.12  is  related  to  the  viscosity  temperature  exponent  for  xenon.  For 
neutral-ion  elastic  interactions,  the  following  cross  sections  measured  by  Miller  et  al.5  are  employed: 

oEL{^e,Xe*  )=  (175.26  -  27.2 log, 0(g))xlO-20m2 


oEL(j(e,Xe**  )=  (103.26  - 17.8  log10(g))x  1(T2IV 

Charge  exchange  concerns  the  transfer  of  one  or  more  electrons  between  an  atom  and  an  ion.  The  cross 
sections  are  assumed  to  follow  the  same  expressions  for  neutral-ion  elastic  collisions.  In  the  present  model,  it  is 
assumed  that  there  is  no  transfer  of  momentum  accompanying  the  transfer  of  the  electron(s).  This  assumption  is 
based  on  the  premise  that  charge  exchange  interactions  are  primarily  at  long  range. 


Plasma  Dynamics 

The  PIC  algorithm  uses  charged  particles  and  determines  the  charge  density  at  the  nodes  of  the  mesh,  based 
on  the  proximity  of  each  particle  to  the  surrounding  nodes.  The  charge  density  is  then  used  to  calculate  the 
potential  at  the  nodes.  The  plasma  potential  can  be  described  using  the  Boltzmann  relationship  that  is  derived  by 
keeping  only  the  dominant  terms  of  the  electron  momentum  equation  and  assuming  isothermal  electrons: 


kT 

<t>-<Pref  =  — ln| 
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The  potential  is  then  differentiated  spatially  to  obtain  the  electric  fields  that  are  used  to  transport  the  ions. 


Boundary  Conditions 

For  the  computations  of  the  Hall  thruster  plume  in  12V,  boundary  conditions  must  be  specified  at  the 
thruster  exit  and  along  all  solid  surfaces  in  the  computational  domain. 

Several  macroscopic  properties  of  the  plasma  exiting  the  thruster  are  required  for  the  computations. 
Specifically,  the  plasma  potential,  the  electron  temperature,  and  for  each  species  the  number  density,  velocity, 
and  temperature  are  required.  These  properties  are  determined  using  an  approach  involving  a  mixture  of 
analysis  and  estimation.  Several  types  of  surfaces  are  included  in  the  computation.  They  are  thruster  walls, 
cryopump  surfaces,  baffles,  and  chamber  walls.  Along  these  walls,  the  plasma  potential  is  set  to  zero.  Any  ions 
colliding  with  the  walls  are  neutralized.  When  particles  hit  the  cryopump  surfaces,  a  fraction  of  the  particles  are 
pumped  away,  which  is  characterized  by  a  sticking  coefficient  (a  value  of  0.8  is  used  in  the  present  study).  For 


particles  scattered  back  into  the  flow  field  from  all  surfaces,  a  diffuse  reflection  is  assumed  characterized  by  the 
surface  temperature. 

Plume  Simulation 

The  12V  chamber  at  AEDC  is  an  Electric  Propulsion  (EP)  test  facility,  whose  height  is  about  12  m.  The 
chamber  has  a  relatively  complex  axi-symmetric  geometry.  As  shown  in  Fig.  1,  the  thruster  is  mounted  on  the 
chamber  axis  and  fired  downward.  There  are  three  baffles  in  the  “waist”  area  and  three  more  near  the  center  of 
the  chamber  bottom.  Two  cryopumps  are  employed  to  pump  the  plume  away.  In  the  present  study,  the  plume  is 
generated  by  a  4.5  kW  class  xenon  Hall  thruster  (the  Aerojet  BPT-4000).  If  the  vacuum  chamber  is  assumed 
empty  when  the  thruster  starts  to  fire,  simulation  (Fig.  2)  shows  that  it  will  take  about  1  second  to  balance  the 
mass  in  the  plume  due 


Figure  1.  Mesh  for  the  12V  chamber.  Figure  2.  Approach  to  steady  state  in  the  simulation 

to  the  thruster  firing  and  the  pumping.  There  is  one  main  reason  for  taking  a  long  time  to  reach  a  steady  state  for 
this  high-speed  plume.  The  cryopump  panels  are  at  very  low  temperature  (20K)  so  that  the  thermal  velocity  of 
particles  reflected  from  the  panels  is  very  small,  which  indicates  that  the  characteristic  speed  of  background 
particles  (reflected  neutrals)  is  about  one  order  of  magnitude  smaller  than  the  neutral  velocity  in  the  near  plume 
and  three  orders  of  magnitude  smaller  than  the  ion  velocity  in  the  nozzle  exit. 

The  relatively  large  speed  at  the  thruster  exit  establishes  a  near  field  plume  in  a  very  short  time  whereas  the 
far  field  plume  requires  a  much  longer  time,  which  can  be  illustrated  by  showing  the  total  number  density  at  two 
different  times  (Fig.  3).  The  ions  establish  themselves  much  more  rapidly  as  indicated  by  the  relatively  minor 
changes  observed  in  Fig.  4  for  these  two  times.  The  more  rapid  convergence  of  the  plasma  component  is  further 
illustrated  in  Fig.  5  that  shows  the  behavior  of  plasma  potential.  The  different  behavior  of  neutrals  and  ions  is 
also  illustrated  in  Fig.  6  by  showing  the  streamlines  of  the  individual  species.  Specifically,  the  ions  are  emitted 
from  the  thruster  exit  and  are  lost  on  any  surfaces,  whereas  the  neutrals  come  from  both  the  thruster  exit  and  the 
surfaces  of  baffle  and  chamber,  and  are  only  removed  by  the  cryopumps. 
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Figure  3.Total  number  density  at  different  times. 
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Figure  4.  Ion  number  density  at  different  times. 
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Plasma  potential  at  different  times. 


left:  neutral  streamlines:  right:  ion  streamlines 
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Figure  6.  Streamlines  of  atoms  and  ions. 


Physics  Modeling 

A  series  of  plume  simulations  are  conducted  to  quantify  the  effects  of  different  physical  mechanisms  on  the 
plasma  plume  structure.  Specifically,  three  simulations  are  performed:  (1)  with  charge  exchange  and  electro¬ 
static  fields  (DSMC,  PIC,  CEX);  (2)  with  electro-static  fields  and  no  charge  exchange  (DSMC,  PIC),  and  (3) 
with  charge  exchange  and  no  electro-static  fields  (DSMC,  CEX).  Results  from  these  simulations  are  compared 
in  Figs.  7a  and  7b  for  the  total  number  density.  In  Fig.  7a,  the  full  simulation  result  is  shown  on  the  left,  and  the 
simulation  omitting  charge  exchange  and  including  electro-static  fields  is  shown  on  the  right.  In  Fig.  7b,  the 
right  hand  solution  includes  charge  exchange  collisions  but  now  omits  the  electric  fields.  There  are  relatively 
minor  differences  between  these  three  solution  results  indicating  that  the  overall  pressure  distribution  in  the  12V 
chamber  is  largely  unaffected  by  the  additional  physics.  The  same  comparisons  are  made  in  Figs.  8a  and  8b  for 
the  plasma  density.  Figure  8a  indicates  that  charge  exchange  has  a  relatively  weak  effect  on  the  plasma 
distribution  except  for  the  backflow  region  behind  the  thruster.  By  contrast,  Fig.  8b  shows  that  omitting  the 
electric  fields  creates  a  significantly  different  plasma  plume  structure 


left:  DSMC'PIC:  right:  DSMC.  no  CEX 


left:  DSMC'PIC:  right:  DSMC,  CEX 
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Figure  7a.  Total  number  density 


obtained  with  different  physical  models. 
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Figure  7b.  Total  number  density 


obtained  with  different  physical  models. 


left:  DSMC  PIC:  right:  DSMC.  no  CEX 


t-  W 
'  ■  '■>  i-.:  % 


o 


0  12  3 

R(m) 

Figure  8a.  Plasma  number  density 
obtained  with  different  physical  models. 
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Figure  8b.  Plasma  number  density 
obtained  with  different  physical  models. 


Comparison  is  also  made  between  the  simulation  results  and  measurement  data.  In  Ref.  7,  the  application  is 
reported  of  a  microwave  interferometry  diagnostic  for  measurement  of  plasma  density  in  12V  in  the  plume  of 
the  BPT-4000  Hall  thruster.  Raw  data  and  best  curve  fits  are  presented  in  Ref.  7  of  the  plasma  density 
distributions.  In  our  use  of  these  data  sets,  the  error  bar  is  set  as  50%  of  the  curve  fitted  value  based  on  our 
observations  of  the  effectiveness  of  the  curve  fit.  In  Fig.  9a,  comparisons  between  the  three  different 
simulations  are  shown  for  the  radial  plasma  density  profiles  at  five  different  axial  distances  from  the  thruster 
exit  plane.  These  profiles  show  quantitatively  the  same  trends  illustrated  in  Figs.  7-8.  Namely,  that  omitting 
the  electric  fields  has  a  much  greater  effect  on  the  plume  structure  than  omitting  the  charge  exchange  collisions. 
In  Fig.  9b,  the  full  DSMC/PIC  simulation  results  (including  both  charge  exchange  and  electric  fields)  are 
compared  with  the  curve-fits  provided  by  Meyer  et  al.7  The  simulation  and  measurement  profiles  are  in 


remarkably  good  agreement  at  all  locations  in  the  plume  considered  The  two  data  sets  lie  within  the  uncertainty 
level  of  the  curve-fit  data  at  all  points.  For  reference,  Fig.  9c  is  taken  directly  from  Ref.  7  and  compares  the 
curve  fits  with  the  raw  measurement  data.  Comparison  of  the  raw  measurement  data  in  Fig.  9c  with  the 
DSMC/PIC  results  in  Fig.  9b  shows  even  better  agreement  than  with  the  curve-fit  data.  The  excellent  levels  of 
agreement  obtained  in  these  comparisons  serves  as  a  strong  validation  of  the  DSMC-PIC  code  developed  for 
analysis  of  plasma  plumes  in  the  12V  facility. 


Figure  9a.  Radial  Profiles  of  plasma  number  density  Figure  9b.  Radial  Profiles  of  plasma  number  density 


Figure  9c.  Radial  Profiles  of  plasma  number  density  (taken  from  Ref.  7). 


Conclusions 

'  A  general  purpose,  hybrid  DSMC-PIC-fluid  code  has  been  developed  for  simulation  of  plasma  plumes  from 
thrusters  operated  in  the  12V  electric  propulsion  facility  at  AEDC.  The  code  was  applied  to  model  the  plasma 
plume  structure  from  the  BPT-4000  Hall  thruster.  A  series  of  simulations  was  performed  in  order  to  assess  the 
sensitivity  of  the  computed  results  to  inclusion  of  different  levels  of  physical  modeling  fidelity.  Comparison  of 
these  results  indicated  that  the  electric  fields  have  a  significantly  stronger  impact  on  the  plasma  plume  structure 
than  charge  exchange  collisions.  The  physical  accuracy  of  the  full  plume  simulation  was  assessed  through 
comparisons  of  previous  measurements  of  plasma  number  density  in  the  thruster  plume  obtained  with  a 
microwave  interferometer.  Excellent  agreement  was  obtained  between  simulation  and  measurement  for  all 
plume  locations  considered.  The  simulation  tool  is  therefore  considered  validated  for  application  to  this  class  of 
Hall  thruster. 
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